#ifndef __CHECK_HPP__
#define __CHECK_HPP__

#include <cmath>  

int fequal(double x, double y) {
  double diff, epsilon = 1e-5;
    int isnan_x, isnan_y; 
    int isinf_x, isinf_y; 

    isnan_x = std::isnan(x);
    isnan_y = std::isnan(y);
    if (isnan_x && isnan_y) return 1;
    if (isnan_x || isnan_y) return 0;

    isinf_x = std::isinf(x);
    isinf_y = std::isinf(y);

    if (isinf_x && isinf_y) {
	return (x > 0 && y > 0) || (x < 0 && y < 0);
    }
    if (isinf_x || isinf_y) return 0;

    diff = fabs(x-y);
    // NTV: 03/30/2009 need this case distinction otherwise comparison against
    // 0.0 will always fail.
    if ((x == 0.0 && y != 0.0) ||
	(x != 0.0 && y == 0.0)) {
	return diff <= epsilon;
    }

    // NTV: 05/06/2009 need this case to rule out differences due to
    // imprecisions close to 0 on cell. Ideally it would only be activated on
    // cell.
    if (x <= epsilon && -epsilon <= x && 
	y <= epsilon && -epsilon <= y) {
	return diff <= epsilon;
    }

    return diff <= epsilon * fabs(x) ||
	diff <= epsilon * fabs(y);
}


int check_matrices(const char * name1, const char * name2,
                   const void * A, const void * B,
                   int m, int n) {
    int i, j;
    const double * X = (const double *)A;
    const double * Y = (const double *)B;
    int errors = 0;
    int max_errors = 10;
    
    for (i = 0; i < m; i++) {
	for (j = 0; j < n; j++) {
	    double x = X[i*n+j];
	    double y = Y[i*n+j];
	    if (! fequal(x, y)) {
		if (++errors >= max_errors) goto EXIT;
		printf("(%s[%d][%d]=%lg) != (%s[%d][%d]=%lg) (diff=%lg)\n",
		       name1, i, j, x,
		       name2, i, j, y,
		       x-y);
	    }
	}
    }
EXIT:

    if (errors > 0) {
	printf("Number of errors=%d (max reported errors %d)\n",
	       errors, max_errors);
    }

    return errors > 0;
}
#endif // __CHECK_HPP__
